Untitled

wd<- "/Users/chengg/Google Drive/EPICON/Mycobiome/Fungal ITS/statistic/Total.fungi"
setwd(wd)

library(vegan)
library(ecodist)

rm(list = ls())
rm(list = ls())

load("EPICON.data.preparation.RC.bNTI.ted.2019.04.19.Rdata")

fung<-fung.rar  #rarefied fungal communities
env<-env        # meta data

plot.gd<-function (habitat, j, k, z) #create a 'plot.gd" function
  #habitat: sample types (Root, Rhizosphere, Soil, Leaf) 
  # k: the y axis value to print mantel result, 0.1 to 0.9
  # the abbreviation of the habitat (Leaf, Root, Rhizo, Soil)
{
  fung1<-fung[env$Habitat==habitat,]
  env1<-env[env$Habitat==habitat,]
  par(mfrow=c(3,6),mar=c(2, 2, 0.5, 0.5))
  for (i in c(j:17))
  {
    fungx<-fung1[ env1$TP== i,]
    envx<-env1[ env1$TP== i,]
    fung.dist<-vegdist(decostand(fungx,"hellinger"),"bray")
    gd<-distance(data.frame(envx$X, envx$Y))
    td<-distance(envx$TP)
    color=rgb(0,0,0,alpha=0.5) 
    plot(jitter(gd),fung.dist,xlab=" ",ylab=" ",,xlim=c(0,80),ylim=c(0,1), col=color, cex=2)
    text(40, 0.9, sprintf("%s-week%s", z, i), col="red", cex =2.5, font =2)
    ma<-mantel(fung.dist~gd, nperm = 999)
    text(40, k, sprintf("R=%.3f, P=%.3f", ma[1], ma[2]), col="blue", cex =2, font =2)
  }
}

plot.gd("Leaf", 1, 0.7, "Leaf")

plot.gd("Root", 1, 0.1, "Root")

plot.gd("Rhizosphere", 1, 0.1, "Rhizo")

plot.gd("Soil", 0, 0.1, "Soil")

rm(list = ls())
rm(list = ls())

load("EPICON.data.preparation.RC.bNTI.ted.2019.04.19.Rdata")

fung<-fung.rar  #rarefied fungal communities
env<-env        # meta data

plot.td<-function (habitat, j, k, z) #create a 'plot.td" function
  #habitat: sample types (Root, Rhizosphere, Soil, Leaf) 
  # k: the y axis value to print mantel result, 0.1 to 0.9
  # the abbreviation of the habitat (Leaf, Root, Rhizo, Soil)
{
  fung1<-fung[env$Habitat==habitat,]
  env1<-env[env$Habitat==habitat,]
  par(mfrow=c(3,6),mar=c(2, 2, 0.5, 0.5))
  for (plot in c(1:6,11:22))
  {
    fungx<-fung1[ env1$plot== plot,]
    envx<-env1[ env1$plot== plot,]
    fung.dist<-vegdist(decostand(fungx,"hellinger"),"bray")
    gd<-distance(data.frame(envx$X, envx$Y))
    td<-distance(envx$TP)
    color=rgb(0,0,0,alpha=0.5) 
    plot(jitter(td),fung.dist,xlab=" ",ylab=" ",xlim=c(0,18),ylim=c(0,1), col=color, cex=2)
    text(9, 0.9, sprintf("%s-plot%s", z, plot), col="red", cex =2.5, font =2)
    ma<-mantel(fung.dist~td, nperm = 999)
    text(9, k, sprintf("R=%.3f, P=%.3f", ma[1], ma[2]), col="blue", cex =2, font =2)
  }
}

plot.td("Leaf", 1, 0.1, "Leaf")

plot.td("Root", 1, 0.1, "Root")

plot.td("Rhizosphere", 1, 0.1, "Rhizo")

plot.td("Soil", 0, 0.1, "Soil")

rm(list = ls())
rm(list = ls())


load("EPICON.data.preparation.RC.bNTI.ted.2019.04.19.Rdata")

fung<-fung.rar[env$TP>0,]  #rarefied fungal communities
env<-env[env$TP>0,]        # meta data
env$Treatment2<-factor(env$Treatment1, labels  = c("CON", "POST", "PRE")) #"CON", "POST", "PRE" are the abbre of Control, Pre_flowering_drought, Post_flowering_drought
par(mfrow=c(3,4),mar=c(2, 2, 0.5, 0.5))
  for (treat in c("CON",  "PRE", "POST")) {
    fung1<-fung[env$Treatment2== treat,]
    env1<-env[env$Treatment2== treat,]
    
    for (habitat in c("Leaf","Root", "Rhizosphere", "Soil")){
      fungx<-fung1[env1$Habitat==habitat ,]
      envx<-env1[env1$Habitat==habitat ,]
      
      fung.dist<-vegdist(decostand(fungx,"hellinger"),"bray")
      gd<-distance(data.frame(envx$X, envx$Y))
      td<-distance(envx$TP)
      color=rgb(0,0,0,alpha=0.5) 
      plot(jitter(td),fung.dist,xlab=" ",ylab=" ",xlim=c(0,18),ylim=c(0,1),  col=color, cex=0.5)
      text(9, 0.95, sprintf("%s-%s", habitat, treat), col="red", cex =2.5, font =2)
      ma<-mantel(fung.dist~td, nperm = 999)
      text(9, 0.05, sprintf("R=%.3f, P=%.3f", ma[1], ma[2]), col="blue", cex =2, font =2)
    }
    }

par(mfrow=c(1,4),mar=c(2, 2, 0.5, 0.5))
for (treat in c("CON")) {
  fung1<-fung[env$Treatment2== treat,]
  env1<-env[env$Treatment2== treat,]
  
  for (habitat in c("Leaf","Root", "Rhizosphere", "Soil")){
    fungx<-fung1[env1$Habitat==habitat ,]
    envx<-env1[env1$Habitat==habitat ,]
    
    fung.dist<-vegdist(decostand(fungx,"hellinger"),"bray")
    gd<-distance(data.frame(envx$X, envx$Y))
    td<-distance(envx$TP)
    color=rgb(0,0,0,alpha=0.5) 
    plot(jitter(td),fung.dist,xlab=" ",ylab=" ",xlim=c(0,18),ylim=c(0,1),  col=color, cex=0.1)
    text(9, 0.95, sprintf("%s", habitat), col="red", cex =2.5, font =2)
    ma<-mantel(fung.dist~td, nperm = 999)
    text(9, 0.05, sprintf("R=%.3f, P=%.3f", ma[1], ma[2]), col="blue", cex =2, font =2)
  }
}